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Properties of three aircraft wake vortex models, Lamb-Oseen, Burnham-Hallock, and 
Proctor are reviewed. These idealized models are often used to initialize the aircraft wake 
vortex pair in large eddy simulations and in wake encounter hazard models, as well as to 
define matched filters for processing lidar observations of aircraft wake vortices. Basic 
parameters for each vortex model, such as peak tangential velocity and circulation strength 
as a function of vortex core radius size, are examined. The models are also compared using 
different vortex characterizations, such as the vorticity magnitude. Results of Euler and 
large eddy simulations are presented. The application of vortex models in the post- 
processing of lidar observations is discussed. 
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Nomenclature 

vortex circulation (m 2 /s) 
initial vortex circulation (m 2 /s) 

5m to 15m circulation (m 2 /s) 

initial vortex pair descent velocity (m/s) 

wing span (m) 

vortex pair separation (m) 

vortex core size (m) 

dimensional Brunt- Vaisala frequency (V 1 ) 
non-dimensional Brunt- Vaisala frequency = Nb^V^ 1 

potential temperature (K) 

temperature ( K) 

eddy dissipation rate (m 2 /s 3 ) 

non-dimensional eddy dissipation rate = (sb 0 / / 3 F 0 _1 

tangential velocity (m/s) 

velocity vector (m/s) 

axial velocity component (m/s) 

lateral/crosswind velocity component (m/s) 

vertical velocity component (m/s) 

pressure (Pa) 

vorticity vector (s' 1 ) 

vorticity components (s' 1 ) 

symmetric (rate-of- strain) tensor 
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Q 

= 

skew- symmetric (spin) tensor 

A 2 

= 

second largest eigenvalue of S_ 2 + Q 2 

CNR 

= 

carrier to noise ratio 

/d 

= 

Doppler shift (Hz) 

k 

= 

optical wavelength ( /urn ) 


= 

lidar system efficiency 

V 

= 

optical frequency (Hz) 

p 

= 

backscatter coefficient ( m x sf x ) 

F h 

= 

excess noise factor 

E 

= 

pulse energy (J) 

A r 

= 

area of receiver aperture ( m 2 ) 

R 

= 

range ( m ) 

c 

= 

speed of light (m/s) 

h 

= 

Planck’s constant (= 6.62xl0" 34 J-s) 

T cg 


coherence time (s) 

Y(t) 

= 

normalized autocorrelation function 

M 

= 

number of independent incoherent averages 


= 

RMS spectral width (Hz) 

FWHM 

= 

full width at half maximum 

SRF 

= 

signal reduction factor 

FOM 

= 

figure of merit 

fs 

= 

sample frequency (MHz) 

A/ 

= 

spectral resolution (MHz) 

</> 

= 

elevation angle (degrees) 


= 

vortex to range gate distance (m) 

AT 

= 

pulse width (s) 


I. Aircraft Vortex Models 

Three idealized vortex models are defined in this section. In the calculations presented below, the domain was 
defined by (y 9 z) e [-50,50] x [0,1 00] m. The vortex was initialized at (y 0 ,z 0 ) = (0,50) m. Unless stated otherwise, 

the aircraft parameters listed in Table 1 were used in this section. The parameters are based on elliptical wing 
loading. 


Table 1: Aircraft Parameters 


Aircraft Parameter 

B747-400 

Wing Span, B (m) 

64.43 

Initial Vortex Pair Separation, bo = Bn/ 4 (m) 

50.60 

Core Radius, r c (m) 

3.75 ~ 0.05821? 

Initial Circulation, r 0 (m 2 /s) 

565 


A. Lamb-Oseen Model 

The Lamb-Oseen vortex is an analytical solution of the Navier-Stokes equations (Appendix A). It has been used 
extensively for initializing large eddy simulations (Hennemann and Holzapfel 2011; Misaka et al. 2012) and in the 
design of lidar matched filters. 
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The tangential velocity, v e as a function of radial distance, r from the vortex center in the Lamb-Oseen model is 
given by (Holzapfel et al. 2000) 


v e (0 = i 1 - ex P (- 1 - 26 (r /r cfl (!) 

where r c is the vortex core radius which is defined as the radius where the tangential velocity is maximum, r 0 is the 
vortex circulation, and r is the distance from vortex center (y 0 ,z 0 ) 


r =^^iy-y of+i^of ■ 

The tangential velocity in Eq. (1) can be converted to Cartesian velocities using the following relations 

v(y,z)= (. sign)v e - 

r 

v ’ 

w(y,z)= - (sign)vg — 


( 2 ) 


(3) 


where sign is used to define the direction of vortex rotation, +1 for clockwise and -1 for counter-clockwise direction. 
The Cartesian components in Eq. (3) can be combined and re-written in terms of v e 


v 6 


Vv 2 +w 2 


(4) 


Figure 1 shows the tangential velocity field for the Lamb-Oseen model. In this example, a B747 was used, and 
the vortex core radius size was set to 3.75m. 
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Figure 1. Lamb-Oseen (LO) Vortex. Tangential Velocity (mis). r c = 3.15m. The resolution of the underlying 
grid is: Aj = Az = 0.5m. B747-400 T 0 = 565m 2 Is. 


B. Burnham-Hallock Model 

The Burnham-Hallock model (Burnham and Hallock 1982) is given by Eq. (5). It is the most widely used model 
for wake vortex applications, which include processing of lidar observations, initialization of large eddy simulations 
(De Visscher et al. 2013), and modeling of aircraft response to wake encounters (Reimer and Vicroy 1996; Schwarz 
et al. 2010). Figure 2 shows the tangential velocity field given by 


Mr) 


2m- r 2 + r} 


( 5 ) 
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Figure 2. Burnham-Hallock (BH) Vortex. Tangential Velocity (mis). r c = 3.75m. The resolution of the 
underlying grid is: Ay = Az = 0.5 m. B747-400 — ' . T 0 = 565m 2 Is. 


C. Proctor Model 

The Proctor vortex model (Proctor 2000) is based on a detailed analysis of lidar observations of vortex tangential 
velocity, and is given by Eq. (6). This model has been used for the initialization of large eddy simulations of wake 
vortex decay in stratified and turbulent environments. Figure 3 shows the tangential velocity field of the Proctor 
model given by 


1.0939- 


v e 


(r)= 


2/zr 


1-expl- 


(- 10(l Ar c / s ) 0 - 75 )] X [l - exp(- 1 252l(r t r c ) 2 )] if r < 1 Ar c 


2 7TT 


1-exp 


(-10(r/5)°- 75 )] 


( 6 ) 


if r>\ Ar c 
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Figure 3. Proctor (FP) Vortex. Tangential Velocity (mis). r c = 3.15m. The resolution of the underlying grid 
is: Ay = Az = 0.5 m. B747-400 — . T 0 = 565m 2 Is. 


The three vortex models described in this paper are not the only vortex models that have been used for wake 
applications; however, they are the most widely used. Other examples include the Winckelmans (Winckelmans et 
al. 2000), the Jacquin, and the Rankine vortex models (Gerz et al. 2002). 
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II. Comparison of Vortex Models 

Basic characteristics of vortex models are compared in this section: 1) magnitude and distribution of tangential 
and vertical velocity components, 2) circulation strength, and 3) tangential velocity magnitude and distribution as a 
function of the vortex core radius size. 

A. Velocity Components 

Comparisons of the tangential and vertical velocity components of the three models along the vortex centerline 
are shown in Figure 4 and Figure 5, respectively. The crosswind velocity, v along the vortex centerline (z = 50m) is 
zero. 



y (m) 


Figure 4. Comparison of the tangential velocity component {mis) along the vortex centerline {z = 50 m\ and is 
the same profile that would be observed through any line through the vortex center - circularly symmetric. 
r c = 3.75 m. B747-400 — — . T 0 = 565 m 2 /s. 
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Figure 5. Comparison of the vertical velocity component {m/s) along the vortex centerline for r c = 3.15m . 
B747-400 — — . T 0 = 565m 1 /s. 


5 

American Institute of Aeronautics and Astronautics 


B. Circulation 

The circulation of a single vortex as a function of radial distance from the vortex center can be defined in terms 
of the tangential velocity 

r(r)= 2nrvQ{r\ (7) 

Eq. (7) for the three vortex models is plotted in Figure 6 using vortex core radius values of 3.15m and 1.28m. 
v e {r) in Eq. (7) was calculated using Eq. (1), (5) and (6) for the three vortex models. The region between 5- 15m 
has been marked in the figure for reference. Average circulation between 5- 15m is often used to characterize wake 
hazard (Burnham and Hallock 2013). Hazard circulation is discussed in Section V. 




Figure 6. Comparison of circulation (m 2 ls) along the vortex centerline. r c = 3.15m (top) and r c = 1.28 m 
(bottom). B747-400 — . r 0 = 565m 2 /s. 

C. Vortex Core Radius 

A wide range of values for the vortex core radius have been suggested in the literature. These values range 
anywhere from 1% of the wing span (Delisi et al. 2003) to 5% of the wing span (Schwarz et al. 2010). Large eddy 
simulations can use values in the range of 7% of the wing span due to computational constraints. The tangential 
velocities for different models and for different values of core radius ranging from 1% to 6.9% of the wing span are 
shown in Figure 7. 
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Figure 7. Tangential velocity {mis) calculation for different core sizes for the three vortex models. Aircraft is 
a B747-400 . T 0 = 565 m 2 /s. B = 64.43/w. 
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D. Vortex Characterizations 


Several methods have been used for vortex characterization (Jiang et al. 2004). The vorticity vector, CD is given 
by 



e x e y e z 

d d d 


dwl dy -dv / dz 


~co x - 

cd = V x V = det 

dx dy dz 

u v w 


du / dz-dw! dx 
dv / dx-du! dy 


CO y 

co z 


where (D x , co y , and co z are the three components of vorticity in x, y and z-direction. In two dimensions, the 

vorticity has only one component (Figure 8 - left panel). The vorticity, VxV is zero around the vortex except near 
the vortex center where a spike in vorticity is observed. The divergence property ( V • V = 0 ) is satisfied by all three 
models. Jeong and Hussain (1995) define the extent of a vortex based on the second largest eigenvalue, X 2 of 
2 2 

S +Q , where S_ is the rate-of- strain tensor and Q is the spin tensor (Appendix B). A vortex is defined as the 

region in which X 2 is negative. Figure 8 (right panel) shows the comparison of X 2 for the three vortex models along 
the centerline. All three models have similar vortex region definition in X 2 . 




Figure 8. Comparison of vorticity, co x (s' 1 ) along the vortex centerline (left panel) and the second largest 
eigenvalue of S 2 + Q 2 , X 2 (right panel) for r c = 5.15m. B747-400. r 0 = 565m 2 /s. 


III. Euler Solutions of Vortex Pair Descent in Stratified Atmosphere 

The descent of an idealized vortex pair in a non-turbulent stratified atmosphere was simulated using a two- 
dimensional Euler code based on a high-resolution wave propagation method (Ahmad and Lindeman 1997). The 
governing equations for atmospheric flows in two dimensions are given by: 


where 


dU dG dH 

1 1 

dt dy dz 


d_Gy | 
dy 


dH y 
dz 


+ iff + S 5 



P 


pv 

U = 

pv 

, G = 

pv 2 +p 


pw 


pvw 


_po_ 


1 

1 


pw 


' 0 ' 

pvw 


0 

2 

, ¥ = 


pw + p 


PS 

wpQ 


_ 0 _ 


( 9 ) 


( 10 ) 


8 

American Institute of Aeronautics and Astronautics 


In Eq. (9)-(10), p is the density of air, v is the velocity component in the y - direction (crosswind direction), w is 
the velocity component in the z-direction, p is the pressure, 6 is the potential temperature, and g is the acceleration 
due to gravity. The system is closed by an equation of state for pressure, 

p=c 0 (p6>y, (ii) 


where C 0 is a constant given by: 


Q) = 


Rr d 

Rd /c v 


( 12 ) 


In the above relations, y is the ratio of specific heats, and R d is the gas constant for dry air. C p and C v are the 
specific heats of air at constant pressure and volume respectively. p 0 is the reference base-state pressure. In these 
simulations, the atmosphere was assumed to be dry, and the source term S was set to zero ( S represents heat sinks 
and sources due to microphysical processes of cloud formation/dissipation and atmospheric radiation transfer). The 
viscous fluxes, G v and H Y were also set to zero, and only the Euler solutions were considered. A high-resolution 
wave propagation method based on flux-wave decomposition (LeVeque 2002; Ahmad and Lindeman 2007) was 
used to compute Eq. (9)-(12). 

The computational domain was defined by (y 9 z) e [-300, 300] x [0,600] m with ts [0,60] 5 -. The mesh had a 

resolution of Ay = Az = 0.3m (2000 x 2000 cells). The wake vortex pair was approximated by the superposition of 
velocity fields of two counter-rotating vortices defined by the idealized vortex models. The vortex pair was 
initialized at ( y 0 ,z 0 ) = (0,300) m. The vortex core radius, r c was set to 2.255 m (-3.5% of the wing span), and the 

initial circulation strength, r 0 was set to 565 m 2 /s (B747-400). The background atmospheric stability was set to zero 
(A* = 0). Open/farfield boundary conditions were used in the lateral and at the top. The bottom boundary was set to 
solid wall. The simulations were run for a final time of 60s. 

Figure 9 shows the computed velocity and pressure fields at time = 60s. The domain maximum and minimum 
velocities and pressure values are compared in Table 2. The Lamb-Oseen model has the highest pressure deficit. 
Also, the maximum and minimum velocity values predicted by the Lamb-Oseen initialization have larger 
magnitudes compared to the Bumham-Hallock and Proctor models. The differences in simulation results using the 
Proctor and the Bumham-Hallock models for initialization were relatively small. The initial descent velocity, V 0 of 
the vortex pair in all three simulations was IJSm/s compared to l.lllm/s estimated from measurements. 


Table 2: Vortex Pair Descent: Comparison of the models after 60s — r c = 2.255 m 


model 

Vmin (Itl/S) 

Vmax (m/s) 

^ ni in (mis) 

^ ni ax (m/s) 

Pmin (Pfl) 

Pm ax (Pd) 

LO 

-13.66 

13.66 

-15.96 

12.32 

-340.73 

12.03 

BH 

-12.01 

12.01 

-14.32 

10.72 

-270.30 

11.97 

FP 

-11.45 

11.45 

-13.76 

10.17 

-248.59 

11.98 
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Figure 9. Vortex pair velocity and pressure fields at time = 60s. v-velocity {mis) is shown in the top panel, w- 
velocity {mis) is shown in the middle panel, and pressure {Pa) is shown in the bottom panel. Aircraft is a 
B747-400 ^ . r 0 = 565 m 2 ls. B = 64.43/w. The vortex core radius, r c = 2.255 m (-3.5% of B). 
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IV. Vortex Pair Descent in Stratified Turbulent Atmosphere 

In this section, the three idealized vortex models were used to initialize large eddy simulations using the 
Terminal Area Simulation System (TASS). TASS (Proctor 1987) is a cloud-resolving model which has been used 
successfully for simulating a wide range of atmospheric flow applications - microburst and wind shear (Proctor and 
Bowles 1992); transport and decay of aircraft wakes in turbulent atmosphere (Proctor 1996; Han et al. 2000a; Han et 
al. 2000b); sensor characterization (Lai et al. 2010), and convection induced turbulence (Proctor et al. 2002; Ahmad 
and Proctor 2011). TASS computes the primitive variables non-hydrostatic equations in three dimensions and is 
capable of resolving flows at multiple spatial (grid resolutions varying from less than \m to 2 km in the horizontal) 
and temporal scales (few seconds in the case of turbulence eddies to hours for long-lived convective phenomena). 
The model solves prognostic equations for potential temperature, water vapor, cloud droplets, ice crystals, rain, 
snow, and hail and includes a microphysics package for cloud and precipitation development. Subgrid scale 
diffusion is parameterized via a Smagorinsky-type turbulence closure (Smagorinsky 1963), and surface layer 
processes are computed based on the Monin-Obukhov similarity theory. The Smagorinsky model implemented in 
TASS includes a rotation correction for vortical flows (Han et al. 2000b) when simulating the decay and transport 
aircraft wakes. 

The TASS model equations are discretized using fourth-order finite-differences in space for the calculation of 
momentum and pressure fields, and the third-order Leonard scheme (Leonard 1995) is used to calculate the transport 
of potential temperature and water vapor. The Klemp-Wilhelmson time- splitting scheme (Klemp and Wilhelmson 
1978) is used for computational efficiency in which the higher- frequency terms given in the left hand side of Eq. 
(13) and (20) are integrated by enforcing the CFL criteria to take into account sound wave propagation due to 
compressibility effects. The remaining terms in Eq. (14), Eq. (20) and Eq. (21) are integrated using a larger time 
step that is appropriate for anelastic and incompressible flows. Non-reflecting Orlanski boundary conditions 
(Orlanski 1976) are imposed at the outflow boundaries. A sixth-order filter is used to damp-out spurious oscillations 
that may arise due to the use of centered-differencing of momentum and pressure terms. The TASS model equations 
with the assumption of a dry atmosphere can be written as follows 


du t + H dp ' 
dt /?o dxf 


duiUj 

dx j 


duj 

+ U 7 - h 


g(H-l)S i3 


1 dry 

+ 

A) dx j 


(13) 


where u x and u } are the velocity components, X{ and Xj are the spatial dimensions, p is the pressure, g is the 
acceleration due to gravity, is the shear stress tensor, and S is the Kronecker delta. H is the density ratio defined 

by 



e_ 


f p ^ Cvlc P 


(14) 


In Eq. (14), C p (= 1004 J K' 1 kg' 1 ) and C v (= 111 J K' 1 kg' 1 ) are the specific heats of air at constant pressure and 
volume, respectively. p 0 is the base-state pressure, p 0 is the base-state density, and O 0 is the base-state potential 
temperature. The horizontal velocities, atmospheric pressure, potential temperature, and density are decomposed 
into initial time- invariant hydrostatic base-states (uq,Vq,Pq,6q,Pq) which are a function of height only and the 
time-variant perturbations ( u\v\p\6\p ') which are a function of both space and time: 


u(x, y, z, t ) = Uq (z) + u'(x, y, z, t ) 

(15) 

v(x, y, z, t) = v 0 (z) + v'(x, y,z,t ) 

(16) 

p(x,y,z,t ) = p 0 (z) + p'(x,y,z,t ) 

(17) 

0{x,y,z,t) = 0 o (z) + 0'(x,y,z,t ) 

(18) 

p(x,y,z,t ) = p 0 (z) + p’(x,y,z,t ) , 

(19) 


The prognostic equation for pressure perturbation is approximated by 
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( 20 ) 


dp' , C P 

dt C v 


P 


du j 



~ Po8 u jdjj ■ 


The conservation of energy is given in terms of potential temperature with the assumption of a dry adiabatic 
atmosphere 


dQ du du j 

— = — + 6 — -. 

dt dxj dxj 


The system is closed by an equation of state for pressure 


p = R dP T • 


( 21 ) 


( 22 ) 


Simulations were run for two different aircraft (B747 and B767). The first step in TASS wake simulations is the 
generation of background homogeneous turbulence. The methodology is described in Han et al. (2000b). A low 
background turbulence field was generated (£* = 0.07) for both simulation sets, and the atmospheric stability was set 
to neutral (A* = 0). The mesh had a resolution of Ax = 2m, and Ay = Az = 1.5 m for both cases. The vortex core 
radius, r c was set to 4.5m in B747 and 3 m in B767 simulations. The choice of r c depends on mesh resolution, and 
the computational requirements for lower values of r c can be restrictive. Periodic boundary conditions were used in 
all directions. For each aircraft, three simulations were run, each initialized with different vortex model. 

Figures (10)-(11) show the T 2 isosurfaces for the two simulation sets at different times. The theoretical value of 
vortex pair time to link, T L based on Crow and Bates (1976) and modified by Sarpkaya et al. (2001) is given by 


T l = 9 


if £*<0.001 


T l = 9.18 - 180£ * 


if 0.001 <£*<0.0121 


Tl 

Tl 


-1.5583 ln(£*) + 0.1556 if 0.0121 < £* < 0.2535 


0.7474 

£ * J 


3/4 


if £* > 0.2535 


(23) 


The time to link for the two cases is compared with the theoretical value in Tables 3-4. Overall the values 
compare well with the theoretical value for all three vortex initializations. The Lamb-0 seen initialization has the 
fastest 71, and the simulations initialized with the Bumham-Hallock model predict values that are closest to the 
theoretical value. A comparison of domain maximum and minimum values of velocity and pressure fields at 60s 
and 120s is given in Tables 5-8. 


Table 3: Time to Link, J L (B747) Table 4: Time to Link, J L (B767) 


Initialization 

T L 

Lamb-Oseen (LO) 

3.82 

Burnham-Hallock (BH) 

4.13 

Proctor (FP) 

3.92 

Theoretical (Eq. 23) 

4.29 


Initialization 

Tl 

Lamb-Oseen (LO) 

3.92 

Burnham-Hallock (BH) 

4.17 

Proctor (FP) 

4.01 

Theoretical (Eq. 23) 

4.29 
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Table 5: Vortex Decay (B747): Comparison of the models after 60s — r c = 4.5 m 


model 

^ m in (m/s) 

^max (m/s) 

^ini n (m/s) 

Wmax (m/s) 

Pmin (Pd) 

Pmax (Pd) 

LO 

-9.363 

9.342 

-11.832 

7.938 

-177.453 

2.691 

BH 

-8.066 

8.359 

-10.334 

6.846 

-133.842 

1.912 

FP 

-8.983 

9.037 

-11.246 

7.511 

-161.874 

2.301 


Table 6: Vortex Decay (B747): Comparison of the models after 120s — r c = 4.5 m 


model 

^ni in (m/s) 

^max (m/s) 

w ni in (m/s) 

Wmax (m/s) 

Pmin (Pd) 

Pmax (Pd) 

LO 

-9.356 

10.424 

-13.866 

8.578 

-152.349 

8.251 

BH 

-7.845 

8.292 

-12.775 

6.726 

-115.079 

6.700 

FP 

-9.934 

9.084 

-13.709 

7.442 

-140.507 

5.826 


Table 7: Vortex Decay (B767): Comparison of the models after 60s - r c = 3m 


model 

V mi „ (m/s) 

Vmax (m/s) 

w ml „ (m/s) 

ww (m/s) 

Pmin (Pd) 

Pmax (Pd) 

LO 

-6.651 

6.759 

-8.994 

5.561 

-87.910 

3.049 

BH 

-6.221 

6.391 

-8.505 

5.183 

-77.393 

2.747 

FP 

-6.547 

6.690 

-8.978 

5.409 

-84.917 

2.864 


Table 8: Vortex Decay (B767): Comparison of the models after 120s -r c = 3m 


model 

^ni in (m/s) 

^max (m/s) 

W ni in (m/s) 

Wmax (m/s) 

Pmin (Pd) 

Pmax (Pd) 

LO 

-7242 

8.416 

-7.890 

6.341 

-70.409 

2.128 

BH 

-5.793 

6.301 

-7.655 

4.731 

-64.298 

1.559 

FP 

-6.647 

6.693 

-7.322 

4.556 

-69.904 

1.904 
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Lamb-Oseen Initialization (120s) 


Lamb-Oseen Initialization (120s) 




Burnham-Hallock Initialization (120s) Burnham-Hallock Initialization (120s) 




Proctor Initialization (120s) Proctor Initialization (120s) 




Figure 10. Top and side views of X 2 isosurface for the three vortex models - Lamb-Oseen (top), Burnham- 
Hallock (middle), and Proctor (bottom). Time = 1205. Ax = 2m , and Ay = Az = 1.5 m. The vortex core radius, 
r c = 4.5 m. B747 ■ — ' . T 0 = 565 m 2 /s. b 0 = 50 m. The vortex pair was initialized at (^q? z o) = (0,500) m. 
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Lamb-Oseen Initialization (60s) 



Lamb-Oseen Initialization (120s) 



Burnham-Hallock Initialization (60s) Burnham-Hallock Initialization (120s) 




Proctor Initialization (60s) 



Proctor Initialization (120s) 
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Figure 11. Top views of A 2 isosurface for the three vortex model initializations - Lamb-Oseen (top), 
Burnham-Hallock (middle), and Proctor (bottom). Time = 60 s in the left column and time = 120s in the right 
column. Ax = 2//z, and Ay = Az = 1 .5///. The vortex core radius, r c = 3. 0m. B767. r 0 = 360/w 2 /s. b {) = 36m. The 
vortex pair was initialized at (jo? z o) = (0,142.5) m. 
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V. Hazard Circulation 


The circulation, T is defined as the line integral of the velocity around a closed loop, C 

r = j) vrf/, (24) 

where v is the fluid velocity. Using Stokes theorem, the line integral in Eq. (24) can be converted into a surface 
integral 

r = £ v • dl = j* (v X v )-dS = [ co ds ■ (25) 

From Eq. (25), the circulation can be interpreted as the flux of vorticity through the surface, S. Figure 12 shows 
the calculation of a>-dS = co x dy dz for the Proctor vortex model. The circulation strengths for different integration 
surface areas are given in Tables 9-10 for two different values of r c . The hazard circulation used in aircraft response 
models is often r 5 . 15m (Burnham and Hallock 2013). Compared to the Bumham-Hallock and Proctor models, the 
values of r 5 .i 5m are much smaller for the Lamb-Oseen model. 


Table 9: Circulation - B747-400 (r 0 = 565 m 2 /s &r c = 3.15m) 


r (m 2 /s) 

Lamb-Oseen (LO) 

Burnham-Hallock (BH) 

Proctor (FP) 

ro-40m 

565.00 

560.07 

564.48 

Co-15m 

565.00 

531.75 

545.20 

r5_i5m 

60.20 

170.20 

113.83 


Table 10: Circulation - B747-400 (T 0 = 565 m 2 /s &r c = 4.5 m) 


T (m 2 /s) 

Lamb-Oseen (LO) 

Burnham-Hallock (BH) 

Proctor (FP) 

Co-40m 

565.00 

557.94 

564.48 

rYi5m 

565.00 

518.34 

545.20 

r5-i5 m 

119.32 

206.23 

143.49 


65 
60 
55 

M 

45 
40 
35 

-40 -20 0 20 40 

y <ni) 

Figure 12. Calculation of co x dydz (/w 2 /s) for the Proctor model. The left panel shows the entire 
computational domain, and the right panel shows a zoomed view of the computational domain for r c = 3.15m. 
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VI. Lidar Measurements 


A. Velocity Measurement 

For a coherent receiver, the precision with which one can estimate velocity in any given range gate depends on 
the signal carrier to noise ratio (CNR), the width of the signal spectrum, and the number of independent 
measurements used to make the estimate. The CNR and the spectral width are influenced by a variety of factors 
including transmit pulse energy and pulse width, aerosol loading and density variation, wind velocity gradients, and 
the signal processing method utilized. The number of independent measurements that can be used for making the 
estimate is typically dictated by sensor requirements. 

Periodogram based signal processing is often used in real time systems for making velocity estimates because of 
computational efficiency. The periodogram is the squared modulus of the discrete Fourier transform of a subset of 
the discrete samples of the return signal from a single lidar pulse. Each subset of samples of the return signal stream 
is referred to as a gate (or range gate). The number of discrete samples used (the length of the gate) is up to the user; 
however, the number of samples is usually chosen based on the range gated window signal coherence time (Eq. 29). 
Velocity precision is improved (i.e., velocity error is reduced) when the number of independent measurements used 
to make the velocity estimate is increased. This is achieved by accumulating periodograms from the same gate or 
adjacent uncorrelated gates (depending on system requirements). Accumulation beats down shot noise fluctuations 
in the signal spectral estimate including fluctuations caused by speckle, and increases the signal detectivity. 

Figure 13 contains two example normalized periodograms: a single realization on the left and the average of 50 
independent realizations on the right. For this data the central bin, Bin 0, represents no wind or a velocity of 0 m/s. 
The spectral resolution of this data is 1.5625 MHz per frequency bin and the optical wavelength is X = 2 pm. A 
velocity estimate is made by converting the measured Doppler shift to velocity using: 

v r = . (26) 

2 

The peak in the spectrum observed at Bin 2 corresponds to a Doppler shift of f D = 3.125 MHz or a velocity of 
-3.125 m/s (winds towards the lidar). To improve accuracy of the estimate, zero padding or a parabolic fit are used 
to better determine the peak location. The actual Doppler shift in this example is 3.5 MHz corresponding to 
fractional Bin 2.24 and a velocity of -3.5 m/s. 
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Figure 13: Normalized single realization periodogram (left) and 50 pulse accumulation (right). Spectral 
resolution of 1.5625 MHz, wavelength of 2 pm and line of sight wind speed of 3.5 m/s towards the lidar. 

A.l CNR 

CNR is defined as the height of the detected signal spectral peak above a unit noise floor and can be expressed as 
(Henderson et al. 2005; Targ et al. 1991): 


cnr =^-— rV— 4— 

F h hu 2 B R 2 


SRF 


( 27 ) 
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where r]' s is the far field system efficiency of the lidar and in the absence of refractive turbulence (in well-designed 
truncated-beam lidar systems r]' s is in the range of 10-15%), F h is the excess noise factor above the shot noise limit 
(in well-designed coherent receivers F h > 0.9), E (J) is the transmitter pulse energy, h = 6.626xl0 -34 (J\s) is 
Planck’s constant, v (Hz) is the optical frequency, T is the one-way path dependent atmospheric transmission, 
/? ( m~ l sf l ) is the path dependent backscatter coefficient, c is the speed of light, B (Hz) is the receiver noise 
equivalent bandwidth, A r (m 2 ) is the area of the receiver aperture, R ( m ) is range, and SRF is the signal reduction 
factor representing inefficiencies associated with diffraction, focus, refractive turbulence, and beam truncation. 

For a Gaussian signal spectrum the “optimal” operating point in terms of the best trade of range resolution and 
CNR is a matched Gaussian integration window, which gives a noise equivalent bandwidth of (Jacob et al. 2009): 

B n =l/r cg , (28) 

where r cg is Goodman’s coherence time for the received signal field defined as (Goodman 1984): 


T cg = \\r(*)f dr , (29) 


where y(r)is the complex autocorrelation function of the received signal. For a Gaussian signal spectrum and a 
matched Gaussian integration window, Goodman’s coherence time is related to the transmit pulse full width at half 
maximum (FWHM), AT, by 


T 


eg 


71 

2 In 2 


AT ~ 1.505 AT 7 • 


(30) 


Of course, one is free to choose any shape and duration for the signal processing integration window. Increasing 
the duration of the window increases CNR up to the intrinsic limit at the cost or range resolution. It is customary, 
because of its simplicity, to use a rectangular integration window. Doing so results in negligible difference in the 
CNR when compared to the matched Gaussian window of the same duration. 

A.2 Velocity Precision 

The minimum velocity error variance that can be achieved with an unbiased optimal estimator can be expressed 
in terms of the narrowband CNR as (Henderson et al. 2005; Van Trees 1968; Doviak and Zrnic 1984): 



(31) 


where 6v (Hz) is the windowed signal spectrum RMS width, M is the number of independent incoherent averages 
(accumulated periodograms) used to make the estimate, and the third term in Eq. (31) is a saturation term that 
accounts for speckle. Assuming a transform limited Gaussian pulse, matched integration window, uniform 
backscatter, and no wind turbulence, the RMS spectral width is related to the pulse FWHM as: 

dV = 2j°I . (32) 

ylln AT 


A.3 Detectivity 

To make a velocity estimate, one must be able to distinguish the peak of the spectrum from the fluctuations in 
the noise background. The detectivity is a measure of the height of signal spectrum peak relative to the fluctuation 
in the background (i.e., noise floor). The detectivity, or figure of merit (FOM), is defined as the ratio of the 
expected signal peak level to the standard deviation of the noise floor (Henderson et al. 2005). As mentioned 
previously, the CNR is the height of the signal peak above the normalized to unity noise floor. The normalized 
noise floor has a standard deviation of 1 /a/M the detectivity is then: 

FOM = y[\ICNR . (33) 
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The FOM speaks to the statistics of the power in any spectral bin and is not a function of any particular system 
parameter. Hence, it is a good metric for defining the maximum working range of a lidar system and useful when 
comparing the performance of one system to another. Empirical results indicate that a F0M~2 results in 
approximately half of the estimates being good estimates and half being anomalies (or false detections due to noise). 
The range corresponding to FOM ~ 2 represents the effective working range of the coherent system (Henderson et 
al. 2005). 

B. Wake Vortex Measurement 

The wake vortex measurement problem differs from the typical wind measurement problem described in the 
previous section in that the aircraft wake velocity flow field contains significant velocity gradients over lengths 
typically shorter than the pulse length. Velocity gradients increase the spectral content of the return signal (i.e., 
increase the signal spectrum RMS width) and reduce the signal spectral peak (i.e., reduce the CNR). By Eq. (31) 
both of these effects work towards reducing the precision of a velocity measurement. Moreover, the velocity 
gradients observed in a wake velocity field typically generate non-Gaussian spectral profiles. Making a velocity 
estimate on wake influenced spectral measurement only estimates the strongest weighted average velocity 
component over the length of the range gate and provides no information of the velocity distribution that produced 
the measured spectra. The addition of an estimate of the spectral width provides some information on the velocity 
distribution, but the actual velocity profile of the wake along a given line of sight is not recoverable from these two 
estimates. 

If one cannot better resolve the velocity profile by using a shorter pulse (or through more advanced signal 
processing techniques) then a better approach is to more effectively utilize the spectral information that is available. 
For example, for a given spectral estimate, instead of only estimating the peak velocity one can also estimate the 
minimum and maximum velocities above some noise threshold. This information provides additional information 
about the shape of the spectrum. This information combined with knowledge of the system and how the signal 
processing influences the spectra (the effects of the pulse width, range gating, pulse averaging, scan rate, spectral 
resolution, etc.,) can then be used to estimate vortex position and circulation strength (Banakh and Smalikho 2013). 

Another approach used for wake vortex estimation is to compare the lidar measured spectra against model 
spectra based on a vortex flow model (Hannon and Thompson 1994). For transverse viewing of aircraft wakes, 
model spectra are constructed based on a vortex tangential velocity model such as those described earlier in this 
paper. The model spectra take into account all the system and signal processing influences that combine to form the 
spectral estimate (Frehlich and Sharman 2005). 

B.l Transverse Viewing Geometry 

The transvers viewing geometry for wake measurements is shown in Figure 14. The lidar is positioned 
transverse to the direction of flight and scans continuously over some portion of the plane perpendicular to the axial 
component of the vortex velocity flow field. As the aircraft moves through the scan, plan aerosols in the atmosphere 
become entrained in the wake. The energy from each transmitted pulse scatters off the aerosols that are positioned 
along the lidar line of sight. The motion of the aerosols Doppler shifts the energy scattered back to the receiver 
affecting the spectral content of the return signal. 

The lidar can only sense that component of motion along its line of sight. For transverse viewing, this means a 
component of the tangential flow field can be sensed. However, what the lidar actually measures will depend on the 
system characteristics and the signal processing and will not simply be equal to the dot product of the tangential 
velocity and the lidar line of sight vector. 

In the results that follow, a lidar system with an optical wavelength of X = 2 pm and a sample frequency 
/ s = 100 MHz, which equates to a velocity bandwidth of Af= ± 25 m/s with real sampling, is used and for the signal 
processing a range gate length of AR ~ 1.5AT, a spectral resolution of A f= 1.5625 MHz and spectral peak detection 
with parabolic fit are used. All signals are ideal (i.e., shot noise and speckle free) for illustrative convenience. 

For simplicity in this discussion, a single vortex with clockwise rotation will be used as shown in Figure 15. The 
vortex will be that of the B747-400 with r 0 = 565 m 2 /s and r c = 3.75 m. It will be centered at a range of 
R = 1.023 km and at an elevation angle of § = 0 degrees. Results will be shown for a single range gate centered at 
1.023 km and for elevations ranging from (|) = -3.01 to +3.01 degrees with line of sight increments of A4> = 0.07 
degrees as depicted in the left panel of Figure 15. These choices conveniently align one line of sight on each side of 
the vortex core with the peak tangential velocity located at one core radius. The elevation span corresponds to 
approximately 54 m along the vortex radial. 
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Figure 14: Transverse viewing geometry. 




Figure 15: Single vortex geometry used in the analysis. Left panel: range gates at fixed range. Right Panel: 
Geometry for determining radial distance to center of range gate. 

From the geometry shown in the right panel of Figure 15 the radius from the center of the vortex to the center of 
each range gate can be shown to be given by: 


= R 


sin(^) 

cos(^/2) 


The model velocity along the lidar line of sight at the center of each range gate is then: 

v r = Vfl(ty)cOs(^/2), 


(34) 


(35) 


where v e is the model tangential velocity. 

B.2 Model Spectra and Measured Velocity Profiles 

The normalized model signal spectra (Frehlich and Sharman 2005) as a function of elevation angle for the B747- 
400 with r 0 = 565 m 2 /s and r c = 3.75 m as a function of elevation are shown in Figure 16 for pulse widths of 
AT=400ns (top row) and A T= 100 ns (bottom row) for the Bumham-Hallock (left), Lamb-Oseen (center) and 
Proctor (right) wake vortex tangential velocity models described earlier. Inspection of the spectra reveal variations 
in the spectral signatures in the vicinity of the core as a function of frequency and elevation when comparing the 
signatures associated with different vortex models (conceptually similar to what is shown in Figure 5). 
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Figure 16: Normalized model signal spectra for B747-400 with r 0 = 565 m 2 /s for pulse widths of AT = 400 ns 
(top row) and AT = 100 ns (bottom row) for the Burnham-Hallock (left), Lamb-Oseen (center), and Proctor 
(right) wake vortex tangential velocity models. 


The spectral peak (fitted) for each elevation is indicated by the solid black curve in each plot, and for 
comparison, the model velocity in frequency space is indicated by the dashed black curve. As can be seen, the 
measured spectral peak and the model spectral peak values differ significantly; however, less so for the shorter pulse 
width. 

The three model velocity profiles (dashed curve) and the resulting lidar measured velocity profiles (solid curve) 
are compared in Figure 17 for pulse widths Ar=400ns (left column) and A T= 100 ns (right column) using the 
spectral peak values from Figure 16 with Eq. (26). As can be seen the lidar measured velocity profiles are highly 
filtered versions of the model profiles and show negligible difference from model to model. 



Figure 17: Comparison of model and measured velocity profiles. Model profiles (solid) and measured profiles 
(dashed) with Burnham-Hallock (blue), Lamb-Oseen (red), and Proctor (black). 

The Burnham-Hallock model and measured velocity profiles for circulation strengths of r 0 = 465 m 2 /s (black), 
r 0 = 565 m 2 /s (blue) and r 0 = 665 m 2 /s are shown in Figure 18 for pulse widths Ar=400ns (left column) and 
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A T= 100 ns (right column). For the longer pulse, there are only minute differences in the measured velocity 
profiles, while for the shorter pulse differences in the measured velocities are more apparent. As one might guess, 
using the measured velocity profile by itself to make an estimate of the circulation strength in noisy conditions is not 
the best approach because of the system response filtering that occurs. 



Figure 18: Burnham-Hallock model (solid) and measured (dashed) velocity profiles for a B747-400 single 
vortex for a pulse width of AT = 400 ns and r 0 = 465 m 2 /s (black), r 0 = 565 m 2 /s (blue) and T 0 = 665 m 2 /s (red). 


A better approach is one based on a more detailed analysis of the spectral content of the return signal. One can 
define velocity envelopes by thresholding the measured signal spectra above some level to estimate minimum and 
maximum velocities to better define the measured velocity distribution for more accurate estimate of the circulation 
strength (Banakh and Smalikho 2013). In another approach, one can compare directly model spectra of different 
circulation strengths with the measured spectral distribution (Hannon and Thompson 1994). Figure 19-20 show the 
model spectra for circulation strengths r 0 = 465 m 2 ls (left), r 0 = 565 m 2 /s (center), and r 0 = 665 m 2 /s (right) for 
pulse widths AT = 400 ns and AT = 1 00 ns respectively. The model velocity profile and the measured velocity 
profile are also shown in frequency space by the black dashed and solid curves respectively. 

Close inspection of the spectra shown below reveal multiple observable variations in the spectral signatures as a 
function of frequency and as a function of elevation when comparing the signatures associated with different 
circulation strengths. Just as multiple independent measurements help improve the precision of a velocity estimate, 
using all of the spectral data will help improve the precision of a wake circulation estimate; however, performance 
will rely on the accuracy of the model velocity flow with actual flows measured in the field. With noisy signals, 
these comparisons are made using techniques that take into consideration the spectral statistics of the signal and 
noise and assess how well the spectral model fits the measured signature (Frehlich and Sharman 2005). 
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Figure 19: Burnham-Hallock Model signal spectra for a B747-400 single vortex for a pulse width of 
AT = 400 ns and r 0 = 465 m 2 /s (left), r 0 = 565 m 2 /s (center) and r 0 = 665 m 2 /s (right). 
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Figure 20: Burnham-Hallock Model signal spectra for a B747-400 single vortex for a pulse width of 
AT = 100 ns and r 0 = 465 m 2 /s (left), r 0 = 565 m 2 /s (center) and r 0 = 665 m 2 /s (right). 


VII. Summary 

In this paper, the properties of three idealized aircraft wake vortex models were described in detail. The effect of 
vortex core radius size on tangential velocities was discussed. The Lamb-Oseen model predicts the highest peak 
tangential velocity and shows the largest variation in peak tangential velocity as a function of vortex core radius. 
The Burnham-Hallock model has the smallest peak tangential velocity, and the Proctor model exhibits the smallest 
variation in peak tangential velocity with change in core radius size. The tangential velocity approaches zero at the 
vortex center and at large distances from the vortex center for all three models. The tangential velocity attains its 
maximum value at the vortex core radius. The vortex circulation is zero at the vortex center for all three models and 
approaches T 0 at r— ► oo. All three models satisfy the divergence property, and the vortex region definition is similar 
in terms of vorticity and ^ 2 - 

The field maximums and minimums in the two-dimensional Euler simulations were larger for Lamb-Oseen 
initializations compared to the other two models. In the three-dimensional TASS simulations, the theoretical value 
of time-to-link, T L was used for comparing the simulations. Faster vortex linking occurred with Lamb-Oseen 
initializations, whereas the T L values predicted by Burnham-Hallock initializations were closest to the theoretical 
value. The hazard circulations, r 5 . 15m were compared for two vortex core sizes. Compared to the Burnham-Hallock 
and Proctor models, the values of r 5 . 15m were much smaller for the Lamb-Oseen model. An aircraft response model 
for wake encounter based on the Lamb-Oseen vortex will underestimate the encounter hazard. 

A brief introduction to the wake measurement problem was presented highlighting the differences between 
velocity estimation and wake estimation. Coherent pulsed lidar is an effective tool for detecting and measuring 
aircraft wake vortices. Wake measurements rely on close scrutiny of the spectral content of the return signal. In 
order to accurately estimate vortex parameters such as circulation strength, position, and core size, a complete 
understanding of how the lidar system parameters and signal processing affect the spectral content of the measured 
wake signature are required. Wake processing of return signals from long pulse lidars are expected to be less 
sensitive to the choice of vortex model used due to a high degree of system response filtering that occurs. Wake 
processing for shorter pulse systems may be more affected by model choice. 


Appendix A: Navier-Stokes Solution of Lamb-Oseen Vortex 

The Lamb-Oseen vortex model is derived in this appendix. The momentum equation is given by 


dV - - 1 - 9- 

+ [V.V IV-F + — Vp-vV 2 V = 0 . 

a v } p 


(A.l) 


Taking the curl of the momentum equation (A.l) 


Vx 


av 

dt 


+ (v.v)v-F + — Vp-vV 2 V = 0 


(A.2) 


From vector algebra, V x Vp = 0 , for a conservative field, Vxf = 0, and, 
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Vx(k-v)k = Vxfiv(K-f)-Fx®)^ vx(f-v)f = (f-v)©-(®-v)f. 
v 2 ) 


Eq. (A.l) simplifies to 


deb 

dt 


■(v.v)a>-(a-v)v = vV 2 a. 


For two-dimensional flows, Eq. (A.3) can be re-written as 


™ + {v.v) 

dt V ’ 


CO = vV CO . 


Ignoring non-linear terms in Eq. (A.4) 


dco _? 

= vV z co . 

dt 


Eq. (A. 5) can be re-written in cylindrical coordinates as 

^2 


dco _ d co v dco 

dt dr 2 r dr 


Using the similarity solution, 


co(r,t)= where 77 = 


Finding the terms in Eq. (A. 6), 


dco 

dco 

drj 

f' 

-V 

-rt 

~di~ 

drj 

dt 

J 

2yfv 

dco 

dco 

dr/ 

f' 

1 

dr 

dr] 

dr 

J 

*Jvt 

d 2 co 


, 1 



dr 2 

= f" 

yjvt 




Substituting Eq. (A. 8) into Eq. (A.6), 


r+l?--- V'=° 

. 2 *1 J 


The solution of Eq. (A.9) is given by: 


/ = ^4 exp 


f 2 ^ 

~ 7 

v 4 , 


+ B . 


Substituting Eq. (A. 10) into Eq. (A.7), 


f 

a>(r,t)=Aex p 


4 vt 


+ B . 


(A.3) 


(A.4) 


(A. 5) 


(A.6) 


(A.7) 


(A. 8) 


(A.9) 


(A. 10) 


(A.l 1) 
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Circulation, T in terms of vorticity is given by 


T = nr (o 


Substitute co(r,t) fromEq. (A. 11) into Eq. (A. 12) 


r (r,t)=A— 2 


nr exp 


V 4vt J 


+ Bnr 


At / = 0, r(r,0)=fi»- 2 =T 0 and at r = 0, r(0,t)= 4 + 5 = 0 

-r 0 


A = -B = ■ 


nr 


Substituting the values of A and B from Eq. (A. 14) into Eq. (A. 13), 

rM=r 0 


f 

f 2^ 

\ 

-exp 

-r 

+ 1 

4 vt 


V 

V J 

) 


L 0 

Eq. (A. 12) can also be written in terms of tangential velocity 

T = 2 TJTVq . 

Substituting T from Eq. (A. 15) into Eq. (A. 16) 

r. 


v e( r >t) = 
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2 nr 


( 

( 2 ^ 

\ 

- exp 

- r 

+ 1 

4vt , 


V 

V ) 

J 


The tangential velocity, v 0 is maximum at r = r c . Differentiating Eq. (A. 17) w.r.tr and setting it to zero, 

r o 


2 nr* 
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~ r c 

2 r c 

_l avn 

~ r c 

4 vt 

i CaU 

4 vt 

4 vt 

V 

V 

J 


v ' ) 
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Simplifying Eq. (A. 18), 


exp 


( 2 V 

r. 


A-'/ 


2 
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. 4w / 


= 1 . 


Solving Eq. (A. 19) for r c 2 / 4vt , 
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4 vt 


rc , L256r . 1 _ 1-256 


4vt 


Substituting — — = from Eq. (A.20) into Eq. (A. 17), 
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(A. 12) 


(A. 13) 


(A. 14) 


(A. 15) 


(A. 16) 


(A. 17) 


(A. 18) 
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(A.20) 


(A.21) 
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Appendix B: The Velocity Gradient Tensor 

The velocity gradient tensor (V ) is given by (see, e.g., Haimes and Kenwright 1999) 
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V can be decomposed into the symmetric tensor ( ¥ ) and the skew-symmetric tensor ( Q ) 


V=S+Q. 


(B.2) 


The symmetric tensor ( ¥ ) is given by 
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(B.3) 


Eq. (B.3) can be written in the following form 
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(B-4) 


Some properties of the symmetric tensor are as follows: 

T 

• Tensor is symmetric (S = S ) 


• The three eigenvalues of S_ are real 

• The trace of the tensor is the divergence of the flow field, and is the same as the trace of the velocity 
gradient tensor 


tr(5)= V-V 


du dv dw 
dx dy dz 


The skew-symmetric tensor ( Q ) is given by 
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